#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "random.h"
#include "rkvector.h"

RkVector_t * RkVector_alloc(int numelem, ...)
{
  RkVector_t * vec;
  int i;
  va_list args;

  vec = (RkVector_t *) malloc(sizeof(*vec));
  if (!vec)
  {
    fprintf(stderr, "fatal: unable to alloc memory for RkVector\n");
    exit(1);
  }
  vec->numelem = numelem;
  
  if (numelem > 0)
    vec->elem = (double *) malloc(sizeof(*vec->elem) * numelem);

  else
  {
    fprintf(stderr, "error: RkVector number of arguments must be > 0\n");
    exit(1);
  }

  if (!vec->elem)
  {
    fprintf(stderr, "fatal: unable to alloc memory for RkVector\n");
    exit(1);
  }

  va_start(args, numelem);

  for (i=0; i<numelem; i++)
    vec->elem[i] = va_arg(args, double);
    
  va_end(args);

  return vec;
}


void RkVector_free(RkVector_t * vec)
{
  if (vec->elem)
    free(vec->elem);

  free(vec);
}

double RkVector_sum(const RkVector_t * vec)
{
  double sum;
  int i;
  
  sum = 0;
  for (i=0; i<vec->numelem; i++)
    sum += vec->elem[i];
  
  return sum;
}

/* assuming sums to 1 */
int RkVector_sample(const RkVector_t * vec)
{
  double desired;
  double cum;
  int i;
  
  desired = RANDOM_DOUBLE;
  cum = 0;
  for (i=0; i<vec->numelem; i++)
  {
    cum += vec->elem[i];
    if (desired < cum)
      return i;
  }
  fprintf(stderr, "error: RkVector doesn't sum to 1, sum = %f\n", cum);
  exit(1);
}

void RkVector_normalize(RkVector_t * vec)
{
  double sum;
  int i;
  
  sum = 0;
  for (i=0; i<vec->numelem; i++)
    sum += vec->elem[i];

  for (i=0; i<vec->numelem; i++)
    vec->elem[i] /= sum;
  
  return;
}

void RkVector_lognormalize(RkVector_t * vec)
{
  double maxval;
  double sum;
  int i;
  
  /* find the maxval */
  maxval = -HUGE_VAL;
  for (i=0; i<vec->numelem; i++)
    if (vec->elem[i] > maxval)
      maxval = vec->elem[i];

    /* if all the elements are -inf we will assign a uniform weight to
     * the elements
     */
    if (-HUGE_VAL == maxval)
    {
      for (i=0; i<vec->numelem; i++)
      {
        vec->elem[i] = 1.0 / ((double) vec->numelem);
      }
    }
    
    /* if there is a +ve infinity in the weights then we will assign
     * uniform weight to all the +ve infinity */
    else if (HUGE_VAL == maxval)
    {
      /* count the number of +inf */
      sum = 0;
      for (i=0; i<vec->numelem; i++)
      {
        if (HUGE_VAL == vec->elem[i])
        {
          vec->elem[i] = 1;
          sum ++;
        }
        else
        {
          vec->elem[i] = 0;
        }
      }
      /* normalize */
      for (i=0; i<vec->numelem; i++)
      {
        vec->elem[i] /= sum;
      }
    }
  
  /* find the sum after subtracting the max and exponentiating */
  else
  {
    sum = 0;
    for (i=0; i<vec->numelem; i++)
    {
      vec->elem[i] = exp(vec->elem[i] - maxval);
      sum += vec->elem[i];
    }

    /* finally, normalize */
    for (i=0; i<vec->numelem; i++)
      vec->elem[i] /= sum;
  }
}

void RkVector_init(RkVector_t * vec, int numelem, double * elem)
{
  vec->numelem = numelem;
  vec->elem = elem;
}

void RkVector_print(const RkVector_t * vec)
{
  int i;

  printf("[ ");
  
  for (i=0; i<vec->numelem; i++)
    printf("%f ", vec->elem[i]);

  printf("]\n");
}
